Coexistence of Self-Organized Criticality and Intermittent Turbulence in the Solar 

Corona 
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An extended data set of extreme ultraviolet images of the solar corona provided by the SOHO 
spacecraft are analyzed using statistical methods common to studies of self-organized criticality 
(SOC) and intermittent turbulence (IT). The data exhibits simultaneous hallmarks of both regimes, 
namely power law avalanche statistics as well as multiscaling of structure functions for spatial 
activity. This implies that both SOC and IT may be manifestations of a single complex dynamical 
process entangling avalanches of magnetic energy dissipation with turbulent particle flows. 

PACS numbers: 05.65+b, 52.35.Ra, 96.60.qe 



SOC and IT represent two paths to dynamical com- 
plexity in driven, extended nonlinear systems. In clas- 
sical fluid turbulence, scaling is often associated with a 
hierarchical structure of eddies extending over the inertial 
range 1], while in SOC, avalanches of localized instabili- 
ties organize the system toward a steady state exhibiting 
long-range correlations up to the system size Q. Some 
authors have argued that SOC and IT are distinct, un- 
related phenomena 0, 0, 0] .while others have suggested 
a fundamental connection p], 0, H, 0, Ho| . In fact, Bak 
and co-workers have speculated that turbulence in the 
limit of high Reynolds number (HRN) may itself be a 
SOC phenomenon [ol Ull [l2| . with dynamical transitions 
observed in HRN turbulence being critical avalanches of 
dissipation. It has also been shown (loj that, in contrast 
to earlier claims [3|, SOC and IT cannot be distinguished 
by analyzing interoccurrence times between bursts: once 
a finite observation threshold (unavoidable in any physi- 
cal measurement) is introduced, even the ordinary BTW 
sandpile exhibits scale free waiting time statistics j 101 ] in- 
dicative of turbulent systems. 

We propose that coexistence of SOC and IT may be 
a generic feature of magnetized astrophysical plasmas. 
One scenario for this to occur is localized instabilities 
HI 14 1 15 1 switching the plasma between frozen and un- 
frozen states. This process resembles stick-slip or depin- 
ning transitions of SOC models 1|| 17 1 and is to some 
extent analogous to rice pile dynamics [lH, where ki- 
netic energy of grains releases stored potential energy. 
Signatures of SOC and IT have also been obtained in 
MHD simulations 0, 0, E| mimicking footpoint mo- 
tions of coronal magnetic loops. However, the explicit 
complementarity between SOC and IT in astrophysical 
observations has not been demonstrated. 

In order to clarify this issue, we present direct observa- 
tional evidence for coexistence of SOC and IT in the mag- 
netized plasma of the solar corona. Using a single high 
resolution data set, we apply two different analysis meth- 



ods - one for analyzing avalanche statistics of the emis- 
sion field and the other for analyzing structure functions 
of the same field. The energy, area and lifetime statistics 
of avalanches detected in this data set obey robust scal- 
ing laws. Unlike previous studies of flare statistics, we 
use a spatiotemporal event detection algorithm compati- 
ble with the usual definition of avalanches in SOC. Next, 
we show that the same data set exhibits multiscaling and 
extended self-similarity (ESS) of higher order structure 
functions - a hallmark of IT phenomena. The observed 
scaling laws show only weak dependence on average so- 
lar activity and were found both at solar minimum (min) 
and maximum (max), indicating that coexistence of SOC 
and IT is a generic characteristic of coronal behavior. 

Dissipation mechanisms in the corona are activated by 
changes in the configuration of its magnetic field. Con- 
vection of magnetic fields leads to radiative transients, 
plasma jets, and explosive events known as flares 22j, 



231 . |24j . The latter are associated with spatially concen- 
trated release of magnetic energy accompanied by local- 
ized plasma heating up to temperatures of 10 7 K, and can 
be observed by short-wavelength light emission. Flares 
tend to appear at irregular times and locations and ex- 
hibit broadband energy, size and lifetime statistics with 
no obvious characteristics scales. This behavior is often 
interpreted as a signature of SOC [H, US EE El 0, El • 
Further, the configuration of magnetic flux in the corona 
has been shown to form a scale-free network with statis- 
tical features that can be captured with a self-organizing 
network model ^30]. Within active regions, the magnetic 
field also has an intermittent spatial structure that re- 
organizes during large flares [311 ]. Our results suggest, 
though, that complementary coexistence of SOC and IT 
is a generic phenomenon not limited to these large events. 

We have studied time series of full-disk digital images 
of the corona taken by the extreme ultraviolet imag- 
ing telescope (EIT) onboard the SOHO spacecraft [24 1 
in the 195 A wavelength band corresponding to the Fe 
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FIG. 1: Two views of the corona's complexity. Upper panel: 
snapshots of high activity coronal regions used to construct 
avalanches (w a = 350). Lower panel: a snapshot of the con- 
tinuous brightness field at the same time instant. 

XII emission at peak coronal temperatures of 1.6 • 10 6 K. 
The data included two observation periods: 06/29/2001 
- 07/28/2001 (3240 images, solar max, average sunspot 
number 64.0) and 10/22/2005- 12/02/2005 (4407 images, 
solar min, average sunspot number 16.3) with a typical 
time resolution of 13.3 min. To reduce optical distortions, 
we studied only central portion of the Sun disk with the 
linear dimensions 1040 x 1040 Mm (256 x 256 pixels with 
5.6 arcsec resolution). The EIT luminosity w(t,r) was 
analyzed as a function of time t and position r on the 
image plane. The dynamics of w(t, r) captures the redis- 
tribution of radiative flux in consecutive EIT frames due 
to a variety of coronal features such as loops and holes, 
mass ejections, etc. For the purpose of this study, we 
treated w(t, r) as a local measure of coronal activity and 
did not filter the data in an attempt to distinguish be- 
tween different types of coronal events. Our analysis was 
based on two alternative approaches allowing a study of 
w(t,r) both as an impulsive avalanching process and a 
continuum turbulent field - as shown in Fig. 1. 

To identify avalanches, we used a spatiotemporal de- 
tection method 13, 33j that resolves concurrent events. 
First, avalanching regions were identified by applying an 
activity threshold w a representing a background EUV 
flux. Contiguous spatial regions with w(r,t) > w a were 
treated as pieces of evolving avalanches. By checking for 
overlap of common pixels between each pair of consec- 
utive EIT frames, we identified a set of 3-dimensional 
spatiotemporal integration domains Aj(z = 1, ..,N) cor- 
responding to each of N individual avalanches. These do- 
mains of contiguous activity in space and time were used 
to compute the lifetimes, Tj = max(i G Aj) — min(< G Aj), 
the radiative emission flux, Ei = J A w(r,t) drdt, as well 
the peak areas, Aj = max(J* A ^ dr), or maximum num- 
ber of pixels in a snapshot of each avalanche. Active 
regions that split after starting at a unique source were 



considered parts of a single avalanche. Active regions 
that merged were considered as separate avalanches, with 
the common "tail" ascribed to the event that started ear- 
lier. Only events that lasted at least two successive time 
frames and were not truncated by the field of view or 
temporal gaps in observations longer than 40 min were 
selected for subsequent analysis. The robustness of the 
obtained statistics was verified by repeatedly running 
the algorithm with substantially different w a . Due to 
large difference in average emission, (w) , during solar min 
and max, w a were defined relative to (w) for each data 
set. Depending on w a , between 4,830 (1,680) and 26,900 
(5,350) coronal events were detected for solar max (min). 

Fig. 2 shows probability distributions for avalanche 
lifetime, total emission flux and peak area for both solar 
max and min. These statistics can be approximately fit- 
ted by the power-law relations p(T)-T" TT , p(E)~E~ TE 
and p(A)^A~ TA with the exponents being almost inde- 
pendent of w a ■ The same exponents have also been found 
using a fluctuating threshold placed at 3 standard devi- 
ations above the average brightness of each image. The 
large-scale rollovers are due to the lack of events whose 
T approaches the maximum available time scale 5.2 • 10 5 s 
given by the ratio between the latitudinal size of the field 
of view and the rotation velocity at the equator. Such 
events tend to cross the field of view and are therefore 
underrepresented in our sample. 
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FIG. 2: Probability distributions of avalanche lifetime, emis- 
sion flux and peak area for solar max and min at four dif- 
ferent thresholds w a = k (w) with k = 0.4 (squares), fc=0.8 
(diamonds), k—1.2 (triangles) and fe=1.6 (circles). The min 
distributions are shifted downward for comparison. 

The exponents observed at solar max (jt = 2.02±0.05, 
t e = 1.70 ± 0.04, t a = 2.53 ± 0.09) and at solar min 
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(r T = 1.96 ± 0.06, t e = 1.66 ± 0.03, t a = 2.51 ± 0.14) 
are indistinguishable within uncertainties. This is also 
true for the exponents z, De and Da defined by the 
relations T ~ l*. E ~ l° E and A ~ 1 Da , where I is 



the linear avalanche scale [47|. The resulting values are 
z = 1.94 ± 0.12, D E = 3.16 ± 0.21 and D A = 1.50 ± 0.10 
for solar max and z = 2.09 ± 0.17, D E = 3.34 ± 0.25 
and Da = 1-52 ± 0.11 for solar min. All values were 
obtained by averaging over four activity thresholds (w a — 
k(w), where k 6 {0.4,0.8,1.2,1.6}) within fixed ranges 
of scales corresponding to the power-law portions of the 
relations. The reported uncertainties are the standard 
errors from this averaging or from the regression estimate 
at individual thresholds, whichever is larger. Up to these 
errors, the exponents satisfy the probability conservation 
relations z(tt — 1) = De{te — 1) = Da(ta — !)• 

All these results support the hypothesis that the 
corona operates in a SOC state. It is worth emphasiz- 
ing that the exponents reported here have been obtained 
using spatiotemporal definition of avalanches which is 
conceptually much closer to measuring avalanches in nu- 
merical simulations than most of the definitions used in 
previous works on flare statistics - except for a few case 
studies [H, 0] focusing on specific coronal conditions. 
The energy distribution exponent te is smaller than 2 
indicating that the coronal heating is dominated by large 
events as opposed to Parker's scenario of nanoflare heat- 
ing 3J]. This conclusion is in agreement with previous 



estimates based on spatiotemporal detection of coronal 
brightenings [23j, l28j and their spatial detection with sub- 
sequent inte grat ion of the emission fluxes over fixed time 
interval 2^, 26[. The exponent Da matches the fractal 
dimension of nanoflares reported in [26j, whereas tt is 
consistent with previous analyses of threshold-dependent 
inter-occurrence times of x-ray bursts measured over the 
whole Sun [3q . as would be expected if total emission 
were a sum over individual avalanches (36|. Its value as 
well as z « 2 indicate that the corona may operate at a 
mean- field limit 37| . We also note that the avalanche ex- 
ponents reported are consistent with the results of forced 
MHD simulation of turbulent coronal heating [l^] ■ 

As our next step, we have analyzed coronal activity as 
a continuum turbulent field. Such analyses are normally 
accomplished by measuring structure functions [38j | for a 
relevant dynamical variable Sz. Taking Sz proportional 
to differences in the scalar field w(r,t) [39(, we define 
the structure functions of order q to be 



S q (l) = (\w(r,t)-w(r + 8r,t)\i) r 



(1) 



Here 5r is the spatial displacement, I = \Sr\, and averag- 
ing indicated by brackets is performed over all positions 
r. Within the inertial range, S q (l) ~ with C(<?) de- 
fined by the turbulent regime under study. In practice, 
scaling directly in S q (l) is often limited in range. This 
is exactly the issue we faced when analyzing the coronal 
images (Fig. 3). In such situations, the ESS method 



is often used which allows to extend the observed range of 
turbulent scaling making it possible to estimate relative 
exponents [48| . The non-trivial behavior of the moments 
is approximately isotropic within uncertainties (see the 
inset in Fig. 4) and so could not be eliminated by project- 
ing the data onto irreducible representations of a lower 
symmetry group [38 1. 
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FIG. 3: Structure functions of coronal activity. Each function 
is normalized by its maximum value. Inset: ESS scaling of 
the structure functions. The horizontal bar shows the range 
used for estimating exponents. 
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FIG. 4: ESS exponents as compared to the Kolmogorov q/2> 
scaling (K41) and hierarchical models of IT defined in the 
text. Inset: solar maximum exponents obtained for horizon- 
tal (X), vertical (Y), and arbitrary (XY) orientations of the 
displacement vector 8r, as a test for isotropy. Error bars in 
both plots show the discrepancy between horizontal and ver- 
tical ("(6) estimates at solar max. 

The functions S q (l) exhibit ESS over almost the en- 
tire range (Fig. 3, inset). The relative ESS exponent 
C(<?)/C(3) (Fig. 4) shows a clear departure from the Kol- 
mogorov law. Both of these features are typical of IT 
systems. To highlight the origin of this intermittency, 
we have tested several analytical fits encompassed by a 
hierarchical model [U |42| : 

C(q) = (q/g)(l -x) + C(l (1 - x/Cy/s) (2) 

The model contains three tuning parameters defined by 
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the relations Sz ~ I 1 / 9 , t e ~ £ s , and C = 3 — £>, where t e 
is the energy transfer time at smallest inertial scales £ and 
D is the dimension of the dissipative structures. The She- 
Leveque (SL) model [HI with vortex filaments {D = 1) 
is obtained substituting g — 3, x — 2 /3 , and C = 2. The 
Iroshnikov-Kraichnan MHD model [44j assumes g = 4, 
a; = 1/2, and C = 1 with dissipative structures inter- 
preted as current sheets and predicts values of the rela- 
tive exponents which for the q range considered here are 
indistinguishable from the SL model. The combination 
g = 3, x = 2/3, and C = 1 gives the Muller-Biskamp 
(MB) model [42j implying that turbulence occurs in a 
3d MHD system with hydrodynamic scaling and sheet- 
like dissipative structures. As Fig. 4 shows, the latter 
provides the best overall description for our data. 

Our main finding is the simultaneous appearance of 
robust signatures of both SOC and IT in a single time 
series of coronal images - including significantly differ- 
ent phases of the solar cycle. These observations can be 
interpreted using a variety of physical scenarios. One of 
them assumes that turbulence in coronal dissipation is 
passively driven by SOC avalanches of dissipating cur- 
rents which modify the ma gnet ic field and shape its in- 
termittent spatial pattern [31]. A more collaborative 
interaction between SOC and MHD turbulence is also 
possible in which SOC avalanches of reconnecting mag- 
netic loops [25L H3| generate inward and outward plasma 
flows [HI and/or cascades of MHD shock waves [46] work- 
ing as sources of turbulence-driven anomalous resistivity 
regions which may affect the avalanches [l3| 

Among the key open problems is the one of the primary 
physical mechanism of the coronal avalanches. A plau- 



sible model has been proposed by Dahlburg et al. [14 1 
who have shown that secondary instabilities needed to 
support the avalanche can be triggered by the misalign- 
ment between the reconnecting flux tubes. The model 
is based on sheet-like dissipative structures compatible 
with the MB turbulent cascade [42j , and predicts rapidly 
evolving instabilities with the energy dissipation growth 
time of about 200 Alfven transit times across the sheets 
14]. These instabilities, reminiscent of toppling events in 
sandpiles 0, can propagate throughout extended coro- 
nal regions and lead to the multiscale explosive release of 
the coronal energy accompanied with SOC and IT signa- 
tures. Observational verification of this mechanism is a 
task for future research. 

We thank V. Abramenko, A. Klimas and A. Pouquet 
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